% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function uup2 = cgns_grid2grid( grid1,ddc_grid1,ubase1,pbase1, ...
           udofs1,pdofs1,uup1, grid2,ubase2,pbase2,udofs2,pdofs2 )
% uup2 = cgns_grid2grid( grid1,ddc_grid,ubase1,pbase1,udofs1,pdofs1,uup1, ...
%                        grid2,ubase2,pbase2,udofs2,pdofs2 )
%
% Interpolate a CG-NS solution from grid1 to grid2. The interpolation
% is done in fact with an approximate L2 projection, where an L2
% projection is taken on each element and then a simple averaging is
% used to enforce interelement continuity.
%
% The input arguments grid1, udofs1, pdofs1 and uup1 can be arrays in
% order to deal with domain decomposition.

 %---------------------------------------------------------------------
 % Set default data
 xi_default = 1.0/(grid1{1}.d+1); % barycenter
 x_default = 0;
 %---------------------------------------------------------------------

 %---------------------------------------------------------------------
 % Interpolate at the quad. nodes
 igi = 1;
 iei = 1;
 % the assumption here is that the same quad. nodes are used for
 % velocity and pressure
 xg   = zeros(grid2.d,ubase2.m);
 uue  = zeros(ubase1.pk,ubase2.m,grid2.d);
 pe   = zeros(pbase1.pk,ubase2.m);
 xi   = zeros(grid2.d,ubase2.m);
 uPHI = zeros(ubase1.pk,ubase2.m);
 pPHI = zeros(pbase1.pk,ubase2.m);
 uu   = zeros(grid2.d,ubase2.m,grid2.ne);
 p    = zeros(        ubase2.m,grid2.ne);
 e2eg = -ones(grid2.ne,2); % used to speedup serach
 tic
 for ie=1:grid2.ne
   if(mod(ie,floor(grid2.ne/100))==0)
     toc
     disp(['Element ',num2str(ie),' of ',num2str(grid2.ne)]);
     tic
   end
   % Gauss points
   xg = grid2.e(ie).b * ubase2.xig;
   for id=1:grid2.d
     xg(id,:) = xg(id,:) + grid2.e(ie).x0(id);
   end
   % locate the elements
   % check whether we have located a neighbour already
   for i=1:grid2.d+1
     ie_neig = grid2.e(ie).ie(i);
     if(ie_neig>0)
       if(e2eg(ie_neig,1)>0)
         iei = e2eg(ie_neig,1); % element index
         igi = e2eg(ie_neig,2); % grid index
	 break
       end
     end
   end
   for l=1:ubase2.m
     [igi,iei,xii] = locate_point(xg(:,l),grid1,ddc_grid1,[igi,iei],1);
     if(iei<0)
       warning(['Unable to locate point ',num2str(xg(:,l)'), ...
             ' in element ',num2str(ie),'; using default value.']);
       xi(:,l) = xi_default;
       uue(:,l,:) = x_default;
       pe(:,l) = x_default;
       igi = 1; iei = 1;
     else
       for id=1:grid2.d
         idx = grid2.d*(udofs1{igi}.dofs(:,iei)-1) + id;
         uue(:,l,id) = uup1{igi}(idx);
       end
       idx = grid2.d*udofs1{igi}.ndofs + pdofs1{igi}.dofs(:,iei);
       pe(:,l) = uup1{igi}(idx);
       xi(:,l) = xii(2:grid2.d+1);
       e2eg(ie,:) = [iei,igi];
     end
   end
   % interpolate
   for i=1:ubase1.pk
     uPHI(i,:) = ev_pol(ubase1.p_s{i},xi);
   end
   for id=1:grid2.d
     uu(id,:,ie) = sum( uue(:,:,id).*uPHI , 1 );
   end
   for i=1:pbase1.pk
     pPHI(i,:) = ev_pol(pbase1.p_s{i},xi);
   end
   p(:,ie) = sum( pe.*pPHI , 1 );
 end
 %---------------------------------------------------------------------

 %---------------------------------------------------------------------
 % Element L2 projections
 % first we need the reference mass matrix (targer grid)
 for i=1:ubase2.pk
   for j=1:ubase2.pk
     M(i,j) = sum( ubase2.wg .* ubase2.p(i,:) .* ubase2.p(j,:) );
   end
 end
 Mi = inv(M);
 % local projector
 for i=1:ubase2.pk
   wgp(i,:) = ubase2.wg .* ubase2.p(i,:);
 end
 uMiwgp = Mi * wgp;
 clear M Mi wgp
 for i=1:pbase2.pk
   for j=1:pbase2.pk
     M(i,j) = sum( ubase2.wg .* pbase2.p(i,:) .* pbase2.p(j,:) );
   end
 end
 Mi = inv(M);
 % local projector
 for i=1:pbase2.pk
   wgp(i,:) = ubase2.wg .* pbase2.p(i,:);
 end
 pMiwgp = Mi * wgp;
 clear M Mi wgp
 % now the element loop
 uup2 = zeros(2*udofs2.ndofs+pdofs2.ndofs,1);
 nnn = zeros(size(uup2));
 for ie=1:grid2.ne
   for id=1:grid2.d
     idx = grid2.d*(udofs2.dofs(:,ie)-1) + id;
     uup2(idx) = uup2(idx) + uMiwgp * uu(id,:,ie)';
     nnn(idx) = nnn(idx)+1;
   end
   idx = grid2.d*udofs2.ndofs + pdofs2.dofs(:,ie);
   uup2(idx) = uup2(idx) + pMiwgp * p(:,ie);
   nnn(idx) = nnn(idx)+1;
 end
 uup2 = uup2./nnn;
 %---------------------------------------------------------------------

return
